2.2 PPAC 信号处理-II

1. Tracking

利用多个PPAC的x，y信息进行traking的原理如图所示,图中红点代表探测器测量的位置。由于PPAC的探测效率，对于每个入射粒子只有部分PPAC的x或y方向可给出位置信息。假设对于某一个入射粒子，有两个以上(含两个)探测器给出x位置信息 (如图:1Ax,1Bx,3x) 时，入射粒子的x-z径迹可由上述x-z点进行线性拟合得到。y-z径迹的处理方法与x-z径迹的处理方法一致。利用上述线性方程，可内插或外推某一z位置的x和y信息。图中空心黑点代表由线性方程计算得到的位置。

setup

本实验中采用的重建策略为：

1) 若F8PPAC3中有有效的位置而且F8PPAC1和F8PPAC2中至少有一层PPAC有有效的位置。

2) 若F8PPAC3没有有效的位置，则

  • a) 这个事件若被 F8PPAC1和F8PPAC2中的3层或3层以上的PPAC探测。
  • b) 若一个事件只有2层PPAC探测到，要求这两层PPAC分别来自两套PPAC1和PPAC2。

束流径迹-利用所有PPAC的信息

setup

2. PPAC 的位置分辨率

准直束刻度方法

  • 通过准直孔将粒子准直，入射到探测器上。此时测得的位置分辨率$\sigma$与探测器本征位置分辨率$\sigma_0$的关系为 $$ \sigma^2=\sigma^2_0+\sigma^2_{geo} $$ ，其中$\sigma_{geo}$为准直孔引入的分辨率（图a）。

  • 上述方法一般用放射源进行，而PPAC探测器的位置分辨与入射粒子的种类以及能量相关，因此由放射源得到的分辨率不能直接用在束流条件下。

  • 束流条件下，如果在探测器前放置准直孔，则由于束流在准直孔上产生散射(图b)，导致$\sigma_{geo}$，不仅与孔的几何条件有关，也与束流的条件相关。

多探测器径迹测量方法(束流)：

  • 选择被全部PPAC(5层)测到的事件进行径迹重建，得到每个探测器的残差$\Delta x^i=x^i-x^i_{track}$的分布，其中$x^i$和$x^i_{track}$分别为测量和拟合得到的位置。径迹外推到靶所在的位置，得到靶点坐标$tx_{track}$。
  • 每个探测器的残差分布的宽度$\sigma^i_{\Delta x}$是所有探测器本征分辨率的函数，即$\sigma^i_{\Delta x}=\sigma^i_{\Delta x}(\sigma^0_0,\sigma^1_0...,\sigma^4_0)$。
  • 探测器的本征分辨率可通过Monte Carlo模拟得到。具体做法如下：对于一组给定分辨率$(\sigma^0_0,\sigma^1_0...,\sigma^4_0)$，每个探测器的位置$x'^{i}$ 由高斯抽样 $Gauss(x^i_{track},\sigma^i_0)$给出。拟合径迹得到每个探测器的拟合位置$x'^i_{track}$,以及残差$\Delta x'^i=x'^i-x'^i_{track}$分布，径迹外推到靶点的坐标为$t_{x'}$。调节不同的分辨率,使得模拟和实验的残差分布符合。这样就能间接算出每一层PPAC的位置分辨率。此时靶点位置残差$\Delta t_{x'}=t_{x'}-t_{x'track}$分布的宽度即为束流在靶上的位置分辨率。

3. PPAC 探测效率

  • 假设有$N$个粒子穿过探测器灵敏面积，探测器实际探测到的数目为$N_{det}$, 则探测效率$\epsilon$为 $$ \epsilon = \frac{N_{det}}{N} $$

  • 在本实验设置下，可以有很多方法求探测器的效率：

    • 例1.ppac的阳极信号的计数作为分母，阴极计数作为分子（图a）。

    • 例2.选择任意两组(i,j)探测器，拟合得到径迹, 记录径迹穿过探测器k灵敏面积的所有事件数目$N_{ij}$。记录在上述条件下探测器k具有正确位置信息的事件数目$N_{kij}$（图b）。 $$ \epsilon = \frac{N_{kij}}{N_{ij}} $$

  • 用上述做法可分别计算探测器的x，y，x-y的效率

示例代码

下面示例代码假设同时有1A，2A，3探测器给出x，y的位置信息。

Branch

$PPACF8[i][j]$ -Calibrated Data(位置已刻度好)

  • 不正确的位置信息用-999表示。
Branch PPAC code
PPACF8[0][0] PPAC 1 Layer A X (mm) xx[0]
PPACF8[0][1] PPAC 1 Layer A Y (mm) yy[0]
PPACF8[0][2] PPAC 1 Layer A Z for X-plan (mm) xz[0]
PPACF8[0][3] PPAC 1 Layer A Z for Y-plan (mm) yz[0]
PPACF8[0][4] PPAC 1 Layer A Anode time (ns)
PPACF8[1][0] PPAC 1 Layer B X (mm)
PPACF8[1][1] PPAC 1 Layer B Y (mm)
PPACF8[1][2] PPAC 1 Layer B Z for X-plan (mm)
PPACF8[1][3] PPAC 1 Layer B Z for Y-plan (mm)
PPACF8[1][4] PPAC 1 Layer B Anode time (ns)
PPACF8[2][0-4] PPAC 2 Layer A * xx[1],yy[1],xz[1],yz[1]
PPACF8[3][0-4] PPAC 2 Layer B * xx2b,yy2b,xz2b,yz2b
PPACF8[4][0-4] PPAC 3 * xx[2],yy[2],xz[2],yz[2]

利用makeclass生成 tracking.h, tracking.C

root -l f8ppac001.root
>tree->MakeClass("tracking");
>.q
编辑 tracking.h和tracking.C 保存
使用方法：
root -l
> .L tracking.C
> tracking t
>t.Loop()
>.q

root -l tracking.root //分析

tracking.h

在traking.h 代码中增加了用户变量，成员函数SetBranch()，TrackInit() 以及SetTrace()定义

//////////////////////////////////////////////////////////
// This class has been automatically generated on
// Wed Mar 11 09:58:39 2020 by ROOT version 6.18/04
// from TTree tree/tree
// found on file: f8ppac001.root
//////////////////////////////////////////////////////////

#ifndef tracking_h
#define tracking_h

#include <TROOT.h>
#include <TChain.h>
#include <TFile.h>

// Header file for the classes stored in the TTree if any.

class tracking {
public :
   TTree          *fChain;   //!pointer to the analyzed TTree or TChain
   Int_t           fCurrent; //!current Tree number in a TChain

// Fixed size dimensions of array or collections stored in the TTree if any.

   // Declaration of leaf types
   Float_t         PPACF8[5][5];
   Float_t         F8PPACRawData[5][5];
   Int_t           beamTrig;
   Int_t           must2Trig;
   Float_t         targetX,targetY;

     //by user
   Double_t xx[3],xz[3],yy[3],yz[3];//1A,2A,3
   Double_t xx2b[2],yy2b[2],xz2b,yz2b;//2B x,y, 0-measured, 1- fitted.
   Double_t dx[3],dy[3];//residual
   Double_t tx,ty;//target position
   Double_t c2nx,c2ny;//chi2/ndf for xfit,yfit

   // List of branches
   TBranch        *b_PPACF8;   //!
   TBranch        *b_F8PPACRawData;   //!
   TBranch        *b_beamTrig;   //!
   TBranch        *b_must2Trig;   //!
   TBranch        *b_targetX;   //!   
   TBranch        *b_targetY;   //!

   tracking(TTree *tree=0);
   virtual ~tracking();
   virtual Int_t    Cut(Long64_t entry);
   virtual Int_t    GetEntry(Long64_t entry);
   virtual Long64_t LoadTree(Long64_t entry);
   virtual void     Init(TTree *tree);
   virtual void     Loop();
   virtual void     SetBranch(TTree *tree);//by user
   virtual void     TrackInit();//by user
   virtual void     SetTrace(TH2D *h,Double_t k,Double_t b,Int_t min,Int_t max);//by user
   virtual Bool_t   Notify();
   virtual void     Show(Long64_t entry = -1);
};

#endif

#ifdef tracking_cxx
tracking::tracking(TTree *tree) : fChain(0) 
{
// if parameter tree is not specified (or zero), connect the file
// used to generate this class and read the Tree.
   if (tree == 0) {
      TFile *f = (TFile*)gROOT->GetListOfFiles()->FindObject("f8ppac001.root");
      if (!f || !f->IsOpen()) {
         f = new TFile("f8ppac001.root");
      }
      f->GetObject("tree",tree);

   }
   Init(tree);
}

tracking::~tracking()
{
   if (!fChain) return;
   delete fChain->GetCurrentFile();
}

Int_t tracking::GetEntry(Long64_t entry)
{
// Read contents of entry.
   if (!fChain) return 0;
   return fChain->GetEntry(entry);
}
Long64_t tracking::LoadTree(Long64_t entry)
{
// Set the environment to read one entry
   if (!fChain) return -5;
   Long64_t centry = fChain->LoadTree(entry);
   if (centry < 0) return centry;
   if (fChain->GetTreeNumber() != fCurrent) {
      fCurrent = fChain->GetTreeNumber();
      Notify();
   }
   return centry;
}

void tracking::Init(TTree *tree)
{
   // The Init() function is called when the selector needs to initialize
   // a new tree or chain. Typically here the branch addresses and branch
   // pointers of the tree will be set.
   // It is normally not necessary to make changes to the generated
   // code, but the routine can be extended by the user if needed.
   // Init() will be called many times when running on PROOF
   // (once per file to be processed).

   // Set branch addresses and branch pointers
   if (!tree) return;
   fChain = tree;
   fCurrent = -1;
   fChain->SetMakeClass(1);

   fChain->SetBranchAddress("PPACF8", PPACF8, &b_PPACF8);
   fChain->SetBranchAddress("F8PPACRawData",  F8PPACRawData, &b_F8PPACRawData);
   fChain->SetBranchAddress("beamTrig", &beamTrig, &b_beamTrig);
   fChain->SetBranchAddress("must2Trig", &must2Trig, &b_must2Trig);
   fChain->SetBranchAddress("targetX",&targetX,&b_targetX);
   fChain->SetBranchAddress("targetY",&targetY,&b_targetY);
   Notify();
}

Bool_t tracking::Notify()
{
   // The Notify() function is called when a new file is opened. This
   // can be either for a new TTree in a TChain or when when a new TTree
   // is started when using PROOF. It is normally not necessary to make changes
   // to the generated code, but the routine can be extended by the
   // user if needed. The return value is currently not used.

   return kTRUE;
}

void tracking::Show(Long64_t entry)
{
// Print contents of entry.
// If entry is not specified, print current entry
   if (!fChain) return;
   fChain->Show(entry);
}
Int_t tracking::Cut(Long64_t entry)
{
// This function may be called from Loop.
// returns  1 if entry is accepted.
// returns -1 otherwise.
   return 1;
}
#endif // #ifdef tracking_cxx

tracking.C

在tracking的基础上进行了修改

#define tracking_cxx
#include "tracking.h"
#include <TH2.h>
#include <TStyle.h>
#include <TCanvas.h>
#include <TF1.h>
#include <TFitResult.h>

void tracking::SetBranch(TTree *tree)
{
  //measured pos
  tree->Branch("xx",&xx,"xx[3]/D");//1A,2A,3
  tree->Branch("xz",&xz,"xz[3]/D");
  tree->Branch("yy",&yy,"yy[3]/D");
  tree->Branch("yz",&yz,"yz[3]/D");

  //difference between measured and calculated -for pos resolution.
  tree->Branch("dx",&dx,"dx[3]/D");
  tree->Branch("dy",&dy,"dy[3]/D");

  //2B x,y
  tree->Branch("xx2b",&xx2b,"xx2b[2]/D");
  tree->Branch("yy2b",&yy2b,"yy2b[2]/D");

  //target x-y
  tree->Branch("tx",&tx,"tx/D");
  tree->Branch("ty",&ty,"ty/D");

  //ch2/ndf for linear fitting.
  tree->Branch("c2nx",&c2nx,"c2nx/D");
  tree->Branch("c2ny",&c2ny,"c2ny/D");

  tree->Branch("beamTrig",&beamTrig,"beamTrig/I");
  tree->Branch("must2Trig",&must2Trig,"must2Trig/I");

  tree->Branch("targetX",&targetX,"targetX");
  tree->Branch("targetY",&targetY,"targetY");  
}

void tracking::TrackInit()
{
  tx=-999;
  ty=-999;

  //1A
  xx[0]=PPACF8[0][0];
  yy[0]=PPACF8[0][1];
  xz[0]=PPACF8[0][2];
  yz[0]=PPACF8[0][3];

  //2A
  xx[1]=PPACF8[2][0];
  yy[1]=PPACF8[2][1];
  xz[1]=PPACF8[2][2];
  yz[1]=PPACF8[2][3];

  //3
  xx[2]=PPACF8[4][0];
  yy[2]=PPACF8[4][1];
  xz[2]=PPACF8[4][2];
  yz[2]=PPACF8[4][3];

  //2B
  xx2b[0]=PPACF8[3][0];
  yy2b[0]=PPACF8[3][1];
  xz2b=PPACF8[3][2];
  yz2b=PPACF8[3][3];

  xx2b[1]=-1000;
  yy2b[1]=-1000;

}

void tracking::SetTrace(TH2D *h,Double_t k,Double_t b,Int_t min,Int_t max){
    if(h==0) return;
    if(min>=max) return;

    for(int i=min;i<max;i++){
        h->Fill(i,(Int_t)(i*k+b));
    }
}

void tracking::Loop()
{
   TH2D *htf8xz=new TH2D("htf8xz","xz trace by ppac",2200,-2000,200,300,-150,150);
   TH2D* htf8yz=new TH2D("htf8yz","yz trace by ppac",2200,-2000,200,300,-150,150);

  TFile *opf=new TFile("tracking.root","recreate");
  TTree *tree=new TTree("tree","ppac traking");
  SetBranch(tree);

   if (fChain == 0) return;
   Long64_t nentries = fChain->GetEntriesFast();
   Long64_t nbytes = 0, nb = 0;
   for (Long64_t jentry=0; jentry<nentries;jentry++) {
      Long64_t ientry = LoadTree(jentry);
      if (ientry < 0) break;
      nb = fChain->GetEntry(jentry);   nbytes += nb;
      TrackInit();
      bool b1a=abs(xx[0])<150 && abs(yy[0])<150;
      bool b2a=abs(xx[1])<150 && abs(yy[1])<150;
      bool b3=abs(xx[2])<150 && abs(yy[2])<150;
      if(!b1a || !b2a || !b3) continue;

      //fit x-z trajectory
      TFitResultPtr r;
      TGraph *grx=new TGraph(3,xz,xx);
      TF1 *fx=new TF1("fx","pol1",-2000,0);
      r=grx->Fit(fx,"SQ");
      xx2b[1]=fx->Eval(xz2b);
      tx=fx->Eval(0);
      SetTrace(htf8xz,fx->GetParameter(1),fx->GetParameter(0),-1800,0);
      for(int i=0;i<3;i++) dx[i]=xx[i]-fx->Eval(xz[i]);
      c2nx=r->Chi2()/r->Ndf();
      delete grx;
      delete fx;

      //fit y-z trajectory      
      TGraph *gry=new TGraph(3,yz,yy);
      TF1 *fy=new TF1("fy","pol1",-2000,0);
      r=gry->Fit(fy,"SQ");
      yy2b[1]=fy->Eval(yz2b);
      ty=fy->Eval(0);
      SetTrace(htf8yz,fy->GetParameter(1),fy->GetParameter(0),-1800,0);
      for(int i=0;i<3;i++) dy[i]=yy[i]-fy->Eval(yz[i]);
      c2ny=r->Chi2()/r->Ndf(); //对于任何程序的自动拟合，原则上都要输出拟合误差进行评估
      delete gry;
      delete fy;

      tree->Fill();
      if(jentry%10000==0) cout<<"processing "<<jentry<<endl;

   }
   htf8xz->Write();
   htf8yz->Write();
   tree->Write();
   opf->Close();
}
In [1]:
%jsroot on
TFile *ipf=new TFile("tracking.root");
TTree *tree=(TTree*) ipf->Get("tree");
TCanvas *c1=new TCanvas("c1","c1");

束流径迹

In [2]:
TH2D *hxz=(TH2D*) ipf->Get("htf8xz");
hxz->Draw("colz");
c1->Draw();
In [3]:
TH2D *hyz=(TH2D*) ipf->Get("htf8yz");
hyz->Draw("colz");
c1->Draw();

束流在靶上投影

  • 假设靶与束流线垂直。
In [4]:
tree->Draw("ty:tx>>htx(120,-60,60,120,-60,60)","must2Trig","colz");
c1->Draw();
In [5]:
tree->Draw("tx:ty>>htx(120,-60,60,120,-60,60)","beamTrig","colz");
c1->Draw();

$\chi^2/Ndf : tx$

残差，chi2/NDF

In [6]:
tree->Draw("dx[0]>>(200,-5,5)");
c1->Draw();//residual
In [7]:
tree->Draw("dx[0]:dx[1]");
c1->Draw();
In [8]:
tree->Draw("c2ny:dy[0]>>hh(40,-10,10,200,0,1000)","","colz");
c1->Draw();//从chi2/ndf图上可看出，部分事件的径迹拟合误差很大，这部分要在后续数据处理中去掉。
In [9]:
tree->Draw("c2ny>>hh(200,0,1000)","","");
c1->Draw();//从chi2/ndf图上可看出，部分事件的径迹拟合误差很大，这部分要在后续数据处理中去掉。
In [10]:
tree->Draw("tx:ty>>(120,-60,60)","c2nx<10 && c2ny<10 && beamTrig ","colz");
c1->Draw();//
In [11]:
tree->Draw("tx:ty>>(120,-60,60,120,-60,60)","(c2nx>20 || c2ny>20) && beamTrig ","colz");
c1->Draw();//

PPAC2B x,y,x-y的探测效率

In [12]:
Long64_t N_track;
Long64_t Nx_det, Ny_det,Nxy_det;
TCut c2btrack="abs(xx2b[1])<100 && abs(yy2b[1])<100";//径迹穿过探测器的灵敏面积
TCut c2bx="abs(xx2b[0])<100";//x面有正确信号
TCut c2by="abs(yy2b[0])<100";//y面有正确信号

1. 计算穿过探测器灵敏面积的粒子数目

In [13]:
tree->Draw("yy2b[1]:xx2b[1]>>(200,-100,100,200,-100,100)",c2btrack,"colz");//fitted
N_track=tree->GetEntries(c2btrack);//得到给定条件下的计数。
cout<<N_track<<endl;
c1->Draw();

232296

2. 计算在上述条件下，探测器的有效信号数目

In [14]:
tree->Draw("yy2b[0]:xx2b[0]>>(200,-100,100,200,-100,100)",c2bx&&c2by&&c2btrack,"colz");
Nx_det=tree->GetEntries(c2bx && c2btrack);
Ny_det=tree->GetEntries(c2by && c2btrack);
Nxy_det=tree->GetEntries(c2bx && c2by && c2btrack);
cout<<Nx_det<<" "<<Ny_det<<" "<<Nxy_det<<endl;
c1->Draw();

216294 215987 201784

3. 计算x, y, x-y 的探测效率

In [ ]:
Double_t ex,ey,exy;
In [16]:
ex=Double_t(Nx_det)/N_track;
ey=Double_t(Ny_det)/N_track;
exy=Double_t(Nxy_det)/N_track;
In [17]:
TString eff;
eff.Form("PPAC2B:\n eff_x=%.2f%%, \n eff_y=%.2f%%,\n eff_x*eff_y=%.2f%%, \n eff_xy=%.2f%%",ex*100,ey*100,ex*ey*100,exy*100);
cout<<eff.Data()<<endl;

PPAC2B:
 eff_x=93.11%, 
 eff_y=92.98%,
 eff_x*eff_y=86.57%, 
 eff_xy=86.87%
In [ ]:
In [ ]: